##################################################################################
# replication code for:
# Polo, Sara and Julian Wucherpfennig. 
# "Trojan Horse, Copycat, or Scapegoat? Unpacking the Refugees-Terrorism Nexus"
# The Journal of Politics
##################################################################################

rm(list=ls(all=TRUE))

setwd("YOURWORKINGDIRECTORY")

library(plm)
library(lmtest)
library(car)
library(stargazer)
library(arm)
library(psych)
library(Formula)
library(lfe)
library(doBy)
library(ggplot2)

### Figure 1
countries <- c("Turkey", "Russia", "India", "Germany", "Hungary", "Italy", "United States", "Poland", 
               "France", "Great Britain", "Australia", "Peru", "Sweden", "South Korea", "Belgium", "Canada", "Japan", "Mexico", "South Africa", "Brazil", "Argentina", "Serbia", "Spain")

y <- c(87, 87, 83, 83, 83, 77, 73, 73, 75, 74, 71, 66, 64, 59, 62, 62, 72, 59, 58, 58, 53, 36, 23)
data <- data.frame(countries, y)

cairo_ps(file ="fg1.eps", width=6*1.618, height=6, onefile = FALSE, fallback_resolution = 600)
par(mar=c(2,6,1,1) + 0.1) 
xx <- barplot(data[order(data[,2],decreasing=FALSE),][,2],names.arg=data[order(data[,2],decreasing=FALSE),][,1], horiz=TRUE, las=1, xlim=c(0,100))
text(y = xx, x = data[order(data[,2],decreasing=FALSE),][,2]+5, label = data[order(data[,2],decreasing=FALSE),][,2], pos = 2, cex = 1, col = "red")
dev.off()



# clustered standard errors
cclust <- function(x) vcovHC(x,  type="sss", cluster="group")

# read country-year data
load("pw_cy.Rdata")

# read dyadic data
load("pw_dy.Rdata")

# formulas
f0 <- Formula( . ~ lngdp + lnpop + pol2 +pol2sq +civilwar + istatewar + extsupport +Wcw +Wtnt )
f2 <- update(f0, lnattack ~ lnref_tnt.l1 + lnref_ntnt.l1 + .)
f2t <- update(f0, lnattack_tnt  ~ lnref_tnt.l1 + lnref_ntnt.l1 + .)
f2p <- update(f0, lnattack_plac  ~ lnref_tnt.l1 + lnref_ntnt.l1+ .)
f2d <- update(f0, lnattack_dom ~ lnref_tnt.l1 + lnref_ntnt.l1+  .)
f2r <- update(f0, lnattack_isref ~ lnref_tnt.l1 + lnref_ntnt.l1 + .)
f2ri <- update(f0, lnattack_isref ~ lnref_tnt.l1 + maxattack5 + lnreftntxmax5 + lnref_ntnt.l1 +lnrefntntxmax5 + .)
f2rv <- update(f0, lnattack_isref ~ lnref_tnt.l1 + lnref_ntnt.l1 + refagainstciv + .)
f2rvi <- update(f0, lnattack_isref ~ lnref_tnt.l1 + lnref_ntnt.l1 + refagainstciv + lnref_tnt.l1:refagainstciv + lnref_ntnt.l1:refagainstciv + .)

# instrumental variable formulas
f2iv2 <- update(f0,  lnattack ~ . | COW + year  | (lnref_tnt.l1 | lnref_ntnt.l1 ~ lnref_iv_tnt + lnref_iv_ntnt) | COW)
f2div2 <- update(f0, lnattack_dom ~ . | COW + year  | (lnref_tnt.l1  |lnref_ntnt.l1 ~ lnref_iv_tnt + lnref_iv_ntnt) | COW)
f2riv2 <- update(f0, lnattack_isref ~ . | COW + year  | ( lnref_tnt.l1 | lnref_ntnt.l1 ~ lnref_iv_tnt + lnref_iv_ntnt ) | COW)
f2fnr <- update(f0, lnattack_fnotref ~ lnref_tnt.l1 + lnref_ntnt.l1 + .)

# coviariate labels for tables
covlab <- c("ln Refugees",
            "ln Refugees w/ TNT ($\\beta_{1}$)", 
            "ln Refugees w/o TNT ($\\beta_{2}$)", 
            "ln GDP p.c.", 
            "ln Population", 
            "Polity2", 
            "Polity2$^{2}$", 
            "Civil war", 
            "Interstate war", 
            "Intervention", 
            "Neighboring civil war", 
            "Neighboring TNT"
)
k <- length(covlab)

# model specifications
model <- "within"
effect <- "twoways"


###############
### Table 1 ###
###############

#### Transnational attacks #####
m2 <- plm(f2, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2 <- coeftest(m2, vcov.= cclust)[,2] 
W2 <- linearHypothesis(m2, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3 <- plm(f2, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3 <- coeftest(m3, vcov.=cclust)[,2] 
W3 <- linearHypothesis(m3, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4 <- plm(f2, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4 <- coeftest(m4, vcov.=cclust)[,2] 
W4 <- linearHypothesis(m4, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

Wald <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", round(c(NA, W2, W3, W4), digits = 3) ),c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
             c("Year fixed effects", "Yes", "Yes", "Yes", "Yes"))

#### Domestic Attacks ####
m2d <- plm(f2d, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2d <- coeftest(m2d, vcov.=cclust)[,2]
W2d <- linearHypothesis(m2d, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3d <- plm(f2d, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3d <- coeftest(m3d, vcov.=cclust)[,2]
W3d <- linearHypothesis(m3d, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4d <- plm(f2d, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4d <- coeftest(m4d, vcov.=cclust)[,2]
W4d <- linearHypothesis(m4d, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

Waldd <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", round(c(NA, W2d, W3d, W4d), digits = 3) ),c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
              c("Year fixed effects", "Yes", "Yes", "Yes", "Yes"))

#### Attacks against against refugees ####
m2r <- plm(f2r, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2r <- coeftest(m2r, vcov.=cclust)[,2]
W2r <- linearHypothesis(m2r, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3r <- plm(f2r, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3r <- coeftest(m3r, vcov.=cclust)[,2]
W3r <- linearHypothesis(m3r, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4r <- plm(f2r, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4r <- coeftest(m4r, vcov.=cclust)[,2]
W4r <- linearHypothesis(m4r, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

Waldr <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", round(c(NA, W2r, W3r, W4r), digits = 3) ),c("Country fixed effects", "Yes", "Yes", "Yes", "Yes"),
              c("Year fixed effects", "Yes", "Yes", "Yes", "Yes"))



WaldFE <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2, W3, W4, W2d, W3d, W4d, W2r, W3r, W4r), digits = 3) )), 
               c("Country fixed effects", rep("\\text{Yes}", 9)),
               c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2,m3,m4, m2d,m3d,m4d ,m2r,m3r,m4r,
          se=list(se2,se3,se4, se2d, se3d, se4d, se2r,se3r,se4r),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldFE,
          covariate.labels = covlab[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Fixed effects regression estimates",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE ,
          font.size = "scriptsize",
          out= "table.FE.tex"
) 


###############
### Table 2 ###
###############

# transnational attacks
m2iv <- felm(f2iv2, data= Tscs)
W2iv <- lfe:::waldtest(m2iv, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

m3iv <- felm(f2iv2, data= Tscs, subset=oecd==0)
W3iv <- lfe:::waldtest(m3iv, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

m4iv <- felm(f2iv2, data= Tscs, subset=oecd==1)
W4iv <- lfe:::waldtest(m4iv, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]


F2iv <- paste( round(condfstat(m2iv)[1],0), round(condfstat(m2iv)[2], 0), sep="/")
F3iv <- paste( round(condfstat(m3iv)[1],0), round(condfstat(m3iv)[2], 0), sep="/")
F4iv <- paste( round(condfstat(m4iv)[1],0), round(condfstat(m4iv)[2], 0), sep="/")

# domestic attacks
m2ivd <- felm(f2div2, data= Tscs)
W2ivd <- lfe:::waldtest(m2ivd, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

m3ivd <- felm(f2div2, data= Tscs, subset=oecd==0)
W3ivd <- lfe:::waldtest(m3ivd, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

m4ivd <- felm(f2div2, data= Tscs, subset=oecd==1)
W4ivd <- lfe:::waldtest(m4ivd, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

# attacks against refugees
m2ivr <- felm(f2riv2, data= Tscs)
W2ivr <- lfe:::waldtest(m2ivr, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

m3ivr <- felm(f2riv2, data= Tscs, subset=oecd==0)
W3ivr <- lfe:::waldtest(m3ivr, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

m4ivr <- felm(f2riv2, data= Tscs, subset=oecd==1)
W4ivr <- lfe:::waldtest(m4ivr, ~lnref_tnt.l1 - lnref_ntnt.l1, type="cluster")[1]

WaldIV <- list(c("F-statistic IV", rep(c(F2iv, F3iv, F4iv),3)),
               c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2iv, W3iv, W4iv, W2ivd, W3ivd, W4ivd, W2ivr, W3ivr, W4ivr), digits = 3) )),
               c("Country fixed effects", rep("\\text{Yes}", 9)),
               c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2iv,m3iv,m4iv, m2ivd,m3ivd,m4ivd ,m2ivr,m3ivr,m4ivr,
          omit.stat=c("f", "adj.rsq", "ser"),
          omit="y",
          add.lines=WaldIV,
          covariate.labels = covlab[-1],
          order=c((k-2):k, 1:(k-1)),
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Instrumental variable regression estimates",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          column.separate = rep(1, 9),
          font.size = "scriptsize",
          out= "table.IV.tex"
)

#################################
### Supplementary Information ###
#################################

### Table A1: descriptive statistics
#### Descriptive statistics
f_all <- update(f0, lnattack  ~ lnattack_dom + lnattack_isref+  lnref.l1 + lnref_tnt.l1 + lnref_ntnt.l1 + .   )

var1 <- covlab
var1[2:3] <- c("ln Refugees w/ TNT" , "ln Refugees w/o TNT")
varnames <- c("ln Transnational attacks", "ln Domestic attacks", "ln Attacks against refugees",  var1, "ln Severe prior attacks")

descr <- model.frame(f_all, Tscs)
stargazer(descr,
          covariate.labels = varnames,
          title= "Descriptive statistics",
          out= "table.DESC.tex"
)


### Figure A1: Trends

ydata.ref <- aggregate(cbind(Tscs$refugees_sum, Tscs$refsum_terrhome, Tscs$refsum_nterrhome, Tscs$attack, Tscs$attack_domestic, Tscs$attack_isref), by=list(Tscs$year), function(x){sum(x, na.rm=TRUE)})
ydata.ref$Group.1 <- as.numeric(as.character(ydata.ref$Group.1))
tab <- t(cbind(ydata.ref$V2/1000000, ydata.ref$V3/1000000))
tab[,24] <- c(NA,NA)

cairo_ps(file ="fgA1.eps", width=8, height=4, onefile = FALSE, fallback_resolution = 600)
par(mar=c(3.1,4.1,3.1,4.1))
barplot(tab , breaks=ydata.ref$Group.1, las=1, ylab="Refugees in millions", bty='n', ylim=c(0,20), names.arg = ydata.ref$Group.1, col = c("grey60", "grey90"))
par(new=T)
plot(ydata.ref$Group.1[1:23], ydata.ref$V5[1:23]+1, axes=F, log = "y", xlim=c(range(ydata.ref$Group.1)), ylim=c(1,20000), xlab="", ylab="", type="l",lty=1, main="",lwd=2, col="red")
axis(side=4,  las=1)
lines(ydata.ref$Group.1[25:47], ydata.ref$V5[25:47]+1, lwd=2, col="red")
lines(ydata.ref$Group.1[1:23], ydata.ref$V4[1:23]+1, lwd=2, col="blue")
lines(ydata.ref$Group.1[1:23], ydata.ref$V6[1:23]+1, lwd=2, col="orange")
lines(ydata.ref$Group.1[25:47], ydata.ref$V4[25:47]+1, lwd=2, col="blue")
lines(ydata.ref$Group.1[25:47], ydata.ref$V6[25:47]+1, lwd=2, col="orange")
mtext("Attacks", side = 4, line = 3, labels=c())
legend(x=1968, y= 30000, legend=c("Transnational attacks", "Domestic attacks", "Attacks against refugees"), col=c("blue", "red", "orange"), lwd=2, bty="n")
legend(x=1993, y= 30000, legend=c("Refugees w/ TNT", "Refugees w/o TNT"), fill= c("grey60", "grey90"), bty="n")
dev.off()

### Table A2: Interaction models

m2ri <- plm(f2ri, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2ri <- coeftest(m2ri, vcov.=cclust)[,2]
W2ri <- linearHypothesis(m2ri, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3ri <- plm(f2ri, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3ri <- coeftest(m3ri, vcov.=cclust)[,2]
W3ri <- linearHypothesis(m3ri, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4ri <- plm(f2ri, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4ri <- coeftest(m4ri, vcov.=cclust)[,2]
W4ri <- linearHypothesis(m4ri, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

stargazer(m2ri,m3ri,m4ri,
          se=list(se2ri,se3ri,se4ri),
          omit.stat=c("f", "adj.rsq"),
          #order = c("lnref.l1", "lnref_tnt.l1", "lnref_ntnt.l1", "lngdp" ,"lnpop" ,"pol2","pol2sq","civilwar", "istatewar", "extsupport", "Wcw", "aid", "armsexp" ),
          order = c(1, 4, 6:14, 2, 3, 5),
          omit="y",
          add.lines=list(               c("Country fixed effects", rep("\\text{Yes}", 3)),
                                        c("Year fixed effects",    rep("\\text{Yes}", 3)) ),
          dep.var.labels = "Attacks against refugees",
          star.cutoffs = c(0.05, 0.01, 0.001),
          covariate.labels = c(covlab[-1], "ln Severe prior attacks", "ln Refugees w/ TNT $\\times$ ln Severe prior attacks", "ln Refugees w/o TNT $\\times$ ln Severe prior attacks"),
          type="latex",
          title = "Interaction models of terrorism against refugees",
          column.labels = c("all countries", "non-OECD", "OECD"),
          #column.separate = c(2, 1, 1),
          font.size= "scriptsize",
          notes = "Clustered standard errors in parentheses",
          out= "table.INTERACTIONS.tex"
) 

### Figure 2: Attacks on refugees conditional on history of severe terrorism

set.seed(321)
prange2 <- c(0.08, 0.5, 0.92) # 84 percent to test for difference at p<.05
k <- 5000 # number of draws

#### model 3ri
draws <- mvrnorm(k, m3ri$coefficients, cclust(m3ri))
# sequence of values of modifying variable
zseq <- seq(min(Tscs$maxattack5[Tscs$oecd==0]), max(Tscs$maxattack5[Tscs$oecd==0]), 0.01)
# empty array to store results
beta1 <- beta2 <- array(NA, c(length(zseq), k))
bs1a <- bs1b <- bs2a <- bs2b <- array(NA,c(length(zseq), 3))
# loop
for (i in 1:k){ 
  b1a <- draws[i, "lnreftntxmax5"]
  b1b <- draws[i, "lnref_tnt.l1"]
  b2a <- draws[i, "lnrefntntxmax5"]
  b2b <- draws[i, "lnref_ntnt.l1"]
  beta1[,i] <- b1b + b1a%*%zseq
  beta2[,i] <- b2b + b2a%*%zseq
}
for (j in 1:length(zseq)){
  bs1b[j,] <- quantile(beta1[j,], probs = prange2)
  bs2b[j,] <- quantile(beta2[j,], probs = prange2)
}
# prepare color
blackalpha <- rgb(0,0,0, alpha=.03)
redalpha <- rgb(1,0,0, alpha=.05)
bluealpha <- rgb(0,0,1, alpha=.05)

cairo_ps(file ="fg2.eps", width=2*6, height=6, onefile = FALSE, fallback_resolution = 600)
par(mfrow=c(1,2))
plot(		c(min(zseq), max(Tscs$maxattack5)),
       c(-.05, .15),  
       main = "non-OECD",  
       xlab = "ln Severe prior attacks",  
       ylab = "Linear predictor",          
       bty="n",
       las=1,
       type = "n")                            

pl1b <- cbind(zseq, bs1b)  #prepare data for plotting
pl2b <- cbind(zseq, bs2b)  #prepare data for plotting

polygon( c(pl1b[,1], rev(pl1b[,1])), c(pl1b[,4] , rev(pl1b[,2])), col = bluealpha, border = 0 )	

polygon( c(pl2b[,1], rev(pl2b[,1])), c(pl2b[,4] , rev(pl2b[,2])), col = redalpha, border = 0 )	
lines(zseq, rowMeans(beta1), lwd=2, col="blue", lty=2)
lines(zseq, rowMeans(beta2), lwd=2, col="red", lty=3)
abline(0,0)
legend("topleft" , legend=c("ln Refugees w/ TNT", "ln Refugees w/o TNT") , bty="n", lty=2:3, lwd=2, col=c("blue","red"))
rug(jitter(Tscs$maxattack5[Tscs$oecd==0], factor=5), col="darkgrey")

draws <- mvrnorm(k, m4ri$coefficients, cclust(m4ri))
# sequence of values of modifying variable
zseq <- seq(min(Tscs$maxattack5[Tscs$oecd==1]), max(Tscs$maxattack5[Tscs$oecd==1]), 0.01)
# empty array to store results
beta1 <- beta2 <- array(NA, c(length(zseq), k))
bs1a <- bs1b <- bs2a <- bs2b <- array(NA,c(length(zseq), 3))
# loop
for (i in 1:k){ 
  b1a <- draws[i, "lnreftntxmax5"]
  b1b <- draws[i, "lnref_tnt.l1"]
  b2a <- draws[i, "lnrefntntxmax5"]
  b2b <- draws[i, "lnref_ntnt.l1"]
  beta1[,i] <- b1b + b1a%*%zseq
  beta2[,i] <- b2b + b2a%*%zseq
}
for (j in 1:length(zseq)){
  bs1b[j,] <- quantile(beta1[j,], probs = prange2)
  bs2b[j,] <- quantile(beta2[j,], probs = prange2)
}

plot(		c(min(zseq), max(Tscs$maxattack5)),
       c(-.05, .15),  
       main = "OECD",  
       xlab = "ln Severe prior attacks",  
       ylab = "Linear predictor",          
       bty="n",
       las=1,
       type = "n")                            

pl1a <- cbind(zseq, bs1a)  #prepare data for plotting
pl1b <- cbind(zseq, bs1b)  #prepare data for plotting
pl2a <- cbind(zseq, bs2a)  #prepare data for plotting
pl2b <- cbind(zseq, bs2b)  #prepare data for plotting

polygon( c(pl1b[,1], rev(pl1b[,1])), c(pl1b[,4] , rev(pl1b[,2])), col = bluealpha, border = 0 )	
polygon( c(pl2b[,1], rev(pl2b[,1])), c(pl2b[,4] , rev(pl2b[,2])), col = redalpha, border = 0 )	
lines(zseq, rowMeans(beta1), lwd=2, col="blue", lty=2)
lines(zseq, rowMeans(beta2), lwd=2, col="red", lty=3)
abline(0,0)
legend("topleft" , legend=c("ln Refugees w/ TNT", "ln Refugees w/o TNT") , bty="n", lty=2:3, lwd=2, col=c("blue","red"))
rug(jitter(Tscs$maxattack5[Tscs$oecd==1], factor=5), col="darkgrey")
dev.off()

### Table A3: controlling for violence by refugees

m2rv <- plm(f2rv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2rv <- coeftest(m2rv, vcov.=cclust)[,2]
W2rv <- linearHypothesis(m2rv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3rv <- plm(f2rv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3rv <- coeftest(m3rv, vcov.=cclust)[,2]
W3rv <- linearHypothesis(m3rv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4rv <- plm(f2rv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4rv <- coeftest(m4rv, vcov.=cclust)[,2]
W4rv <- linearHypothesis(m4rv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2rvi <- plm(f2rvi, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2rvi <- coeftest(m2rvi, vcov.=cclust)[,2]
W2rvi <- linearHypothesis(m2rvi, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3rvi <- plm(f2rvi, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3rvi <- coeftest(m3rvi, vcov.=cclust)[,2]
W3rvi <- linearHypothesis(m3rv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4rvi <- plm(f2rvi, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4rvi <- coeftest(m4rvi, vcov.=cclust)[,2]
W4rvi <- linearHypothesis(m4rvi, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]



WaldRV <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", round(c(W2rv, W3rv, W4rv,W2rvi, W3rvi, W4rvi), digits = 3) ),
                 c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                 c("Year fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"))

stargazer(m2rv,m3rv,m4rv,m2rvi,m3rvi,m4rvi,
          se=list(se2rv,se3rv,se4rv,se2rvi,se3rvi,se4rvi),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldRV,
          covariate.labels = c(covlab[-1], "Violence by refugees", "ln Refugees w/ TNT $\\times$ violence by refugees", "ln Refugees w/o TNT $\\times$ violence by refugees"),
          order = c(1, 2 , 4:12, 3, 13, 14),
          dep.var.labels = c("Attacks against refugees", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          font.size = "scriptsize",
          title = "Results controlling for violence by refugees",
          column.labels = c("all countries", "non-OECD", "OECD", "all countries", "non-OECD", "OECD"),
          notes = "Clustered standard errors in parentheses",
          out= "table.RV.tex"
)


### Table A4: placebo outcomes

m2p <- plm(f2p, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2p <- coeftest(m2p, vcov.=cclust)[,2]
W2p <- linearHypothesis(m2p, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3p <- plm(f2p, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3p <- coeftest(m3p, vcov.=cclust)[,2]
W3p <- linearHypothesis(m3p, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4p <- plm(f2p, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4p <- coeftest(m4p, vcov.=cclust)[,2]
W4p <- linearHypothesis(m4p, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2fnr <- plm(f2fnr, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2fnr <- coeftest(m2fnr, vcov.=cclust)[,2]
W2fnr <- linearHypothesis(m2fnr, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3fnr <- plm(f2fnr, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3fnr <- coeftest(m3fnr, vcov.=cclust)[,2]
W3fnr <- linearHypothesis(m3fnr, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4fnr <- plm(f2fnr, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4fnr <- coeftest(m4fnr, vcov.=cclust)[,2]
W4fnr <- linearHypothesis(m4fnr, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


WaldPLAC <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", round(c(W2p, W3p, W4p,W2fnr, W3fnr, W4fnr), digits = 3) ),
                 c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                 c("Year fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"))

stargazer(m2p,m3p,m4p,m2fnr,m3fnr,m4fnr,
          se=list(se2p,se3p,se4p,se2fnr,se3fnr,se4fnr),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldPLAC,
          covariate.labels = covlab[-1],
          dep.var.labels = c("Transnational attacks (placebo)", "Attacks against refugees (placebo)"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          font.size = "scriptsize",
          title = "Placebo outcome regression estimates",
          column.labels = c("all countries", "non-OECD", "OECD", "all countries", "non-OECD", "OECD"),
          notes = "Clustered standard errors in parentheses",
          out= "table.PLAC.tex"
)

### Table A5: Dyadic Zero Stage
f0iv <- Formula( log(refugees_sum+1)~
                   ( log(capdist) +comlang_ethno ) : ( war2not1  +pol2.1) + ( war2not1 +pol2.1) | year+ dyadid )
iv0 <- felm(f0iv, data=ddata1, na.action = na.omit)
stargazer(iv0, 
          omit.stat=c("f",  "ser", "adj.rsq"),
          covariate.labels = c("War (not involving destination)", "Polity2$_{d}$", "War $\\times$ ln distance", "Polity2$_{d}$ $ \\times$ ln distance",  "War $\\times$ common language", "Polity2$_{d}$ $\\times$ common language"),
          dep.var.labels = c("ln Refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "(Dyadic) ``zero stage'' model",
          add.lines = list(   c("Dyad fixed effects", "Yes"),
                              c("Year fixed effects", "Yes")),
          #notes = "Clustered standard errors in parentheses",
          out= "table.IV.ZERO.tex"
)


### Table A6: IV First Stage Models

f2iv1a <- update(f0,  lnref_tnt.l1 ~ . + lnref_iv_tnt + lnref_iv_ntnt | COW + year | 0 |  COW)
f2iv1b <- update(f0,  lnref_ntnt.l1 ~ . + lnref_iv_tnt + lnref_iv_ntnt | COW + year | 0 |  COW)

m2.1a <- felm(f2iv1a, data= Tscs)
m2.1b <- felm(f2iv1b, data= Tscs)
m3.1a <- felm(f2iv1a, data= Tscs, subset=oecd==0)
m3.1b <- felm(f2iv1b, data= Tscs, subset=oecd==0)
m4.1a <- felm(f2iv1a, data= Tscs, subset=oecd==1)
m4.1b <- felm(f2iv1b, data= Tscs, subset=oecd==1)

covlab.iv <- c("ln Refugees w/ TNT (projected)" , "ln Refugees w/o TNT (projected)" , covlab[3:k])
#dvlab.iv <- c("ln Refugees" , "CHANGE THIS MANUALLY AT THE END", rep(c("ln Refugees w/ TNT","ln Refugees w/o TNT"),3))

Fiv <- c(round(condfstat(m2iv)[1], 0), 
         round(condfstat(m2iv)[2], 0),  
         round(condfstat(m3iv)[1], 0), 
         round(condfstat(m3iv)[2], 0), 
         round(condfstat(m4iv)[1], 0), 
         round(condfstat(m4iv)[2], 0))
k <- length(covlab)

FEs2 <- list( c("Conditional F-statistic", Fiv),  
              c("Country fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
              c("Year fixed effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes") )
stargazer(m2.1a, m2.1b, m3.1a, m3.1b, m4.1a, m4.1b,
          omit.stat=c("f",  "ser" , "adj.rsq"),
          order=c((k-2):k, 1:(k-2)),
          covariate.labels = covlab.iv[-1],
          dep.var.labels = c(rep(c("ln R w/ T","ln R w/o T"),3)),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          #type="text",
          title = "IV First Stage Models",
          column.labels = c("all countries", "non-OECD", "OECD"),
          column.separate = c(2, 2, 2),
          add.lines = FEs2,
          notes = "Clustered standard errors in parentheses",
          font.size = "scriptsize",
          out= "table.IV.FIRST.tex"
)

### Table A7: Terrorist fatalities
f2.k <- update(f2, log(nkill+1) ~ .)

m2.k <- plm(f2.k, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.k <- coeftest(m2.k, vcov.= cclust)[,2] 
W2.k <- linearHypothesis(m2.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.k <- plm(f2.k, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.k <- coeftest(m3.k, vcov.=cclust)[,2] 
W3.k <- linearHypothesis(m3.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.k <- plm(f2.k, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.k <- coeftest(m4.k, vcov.=cclust)[,2] 
W4.k <- linearHypothesis(m4.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


f2d.k <- update(f2d, log(nkill_domestic+1) ~ .)

m2d.k <- plm(f2d.k, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2d.k <- coeftest(m2d.k, vcov.=cclust)[,2]
W2d.k <- linearHypothesis(m2d.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3d.k <- plm(f2d.k, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3d.k <- coeftest(m3d.k, vcov.=cclust)[,2]
W3d.k <- linearHypothesis(m3d.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4d.k <- plm(f2d.k, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4d.k <- coeftest(m4d.k, vcov.=cclust)[,2]
W4d.k <- linearHypothesis(m4d.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


f2r.k <- update(f2d, log(nkill_isref +1) ~ .)

m2r.k <- plm(f2r.k, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2r.k <- coeftest(m2r.k, vcov.=cclust)[,2]
W2r.k <- linearHypothesis(m2r.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3r.k <- plm(f2r.k, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3r.k <- coeftest(m3r.k, vcov.=cclust)[,2]
W3r.k <- linearHypothesis(m3r.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4r.k <- plm(f2r.k, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4r.k <- coeftest(m4r.k, vcov.=cclust)[,2]
W4r.k <- linearHypothesis(m4r.k, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]



WaldKILL <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.k, W3.k, W4.k, W2d.k, W3d.k, W4d.k, W2r.k, W3r.k, W4r.k), digits = 3) )), 
                 c("Country fixed effects", rep("\\text{Yes}", 9)),
                 c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.k,m3.k,m4.k, m2d.k,m3d.k,m4d.k ,m2r.k,m3r.k,m4r.k,
          se=list(se2.k,se3.k,se4.k, se2d.k, se3d.k, se4d.k, se2r.k,se3r.k,se4r.k),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldKILL,
          covariate.labels = covlab[-1],
          dep.var.labels = c("Transnational kills", "Domestickills", "Kills of refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Terrorist fatalities",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.KILL.tex"
) 


### Table A8: Imputed Refugess Data
f2.m <- update(f0, lnattack ~ lnref_tnt_m2.l1 + lnref_ntnt_m2.l1 + .)
f2.d.m <- update(f0, lnattack_dom ~ lnref_tnt_m2.l1 + lnref_ntnt_m2.l1 + .)
f2.r.m <- update(f0, lnattack_isref ~ lnref_tnt_m2.l1 + lnref_ntnt_m2.l1 + .)

m2.m <- plm(f2.m, data=subset(Tscs, mar==1), index = c("COW","year"),  model = model, effect=effect)
se2.m <- coeftest(m2.m, vcov.=cclust)[,2] 
W2.m <- linearHypothesis(m2.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.m <- plm(f2.m, data=subset(Tscs, oecd==0 & mar==1), index = c("COW","year"),  model = model, effect=effect)
se3.m <- coeftest(m3.m, vcov.=cclust)[,2] 
W3.m <- linearHypothesis(m3.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.m <- plm(f2.m, data=subset(Tscs, oecd==1 & mar==1), index = c("COW","year"),  model = model, effect=effect)
se4.m <- coeftest(m4.m, vcov.=cclust)[,2] 
W4.m <- linearHypothesis(m4.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.d.m <- plm(f2.d.m, subset(Tscs, mar==1), index = c("COW","year"),  model = model, effect=effect)
se2.d.m <- coeftest(m2.d.m, vcov.=cclust)[,2] 
W2.d.m <- linearHypothesis(m2.d.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.d.m <- plm(f2.d.m, data=subset(Tscs, oecd==0 & mar==1), index = c("COW","year"),  model = model, effect=effect)
se3.d.m <- coeftest(m3.d.m, vcov.=cclust)[,2] 
W3.d.m <- linearHypothesis(m3.d.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.d.m <- plm(f2.d.m, data=subset(Tscs, oecd==1 & mar==1), index = c("COW","year"),  model = model, effect=effect)
se4.d.m <- coeftest(m4.d.m, vcov.=cclust)[,2] 
W4.d.m <- linearHypothesis(m4.d.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.r.m <- plm(f2.r.m, subset(Tscs, mar==1), index = c("COW","year"),  model = model, effect=effect)
se2.r.m <- coeftest(m2.r.m, vcov.=cclust)[,2] 
W2.r.m <- linearHypothesis(m2.r.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.r.m <- plm(f2.r.m, data=subset(Tscs, oecd==0 & mar==1), index = c("COW","year"),  model = model, effect=effect)
se3.r.m <- coeftest(m3.r.m, vcov.=cclust)[,2] 
W3.r.m <- linearHypothesis(m3.r.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.r.m <- plm(f2.r.m, data=subset(Tscs, oecd==1 & mar==1), index = c("COW","year"),  model = model, effect=effect)
se4.r.m <- coeftest(m4.r.m, vcov.=cclust)[,2] 
W4.r.m <- linearHypothesis(m4.r.m, c("lnref_tnt_m2.l1=lnref_ntnt_m2.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

covlab.m <- covlab
covlab.m[2:3] <- c("ln Refugees w/ TNT, int. ($\\beta_{1}$)", "ln Refugees w/o TNT, int. ($\\beta_{2}$)")

WaldMARBACH <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.m, W3.m, W4.m, W2.d.m, W3.d.m, W4.d.m, W2.r.m, W3.r.m, W4.r.m), digits = 3) )), 
                    c("Country fixed effects", rep("\\text{Yes}", 9)),
                    c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.m,m3.m,m4.m, m2.d.m,m3.d.m,m4.d.m ,m2.r.m,m3.r.m,m4.r.m,
          se=list(se2.m,se3.m,se4.m, se2.d.m, se3.d.m, se4.d.m, se2.r.m,se3.r.m,se4.r.m),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldMARBACH,
          covariate.labels = covlab.m[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Interpolated refugees data",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.MARBACH.tex"
) 


### Table A9: Post 9/11

m2.911 <- plm(f2, data=subset(Tscs, p911==1), index = c("COW","year"),  model = model, effect=effect)
se2.911 <- coeftest(m2.911, vcov.=cclust)[,2] 
W2.911 <- linearHypothesis(m2.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.911 <- plm(f2, data=subset(Tscs, oecd==0 & p911==1), index = c("COW","year"),  model = model, effect=effect)
se3.911 <- coeftest(m3.911, vcov.=cclust)[,2] 
W3.911 <- linearHypothesis(m3.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.911 <- plm(f2, data=subset(Tscs, oecd==1 & p911==1), index = c("COW","year"),  model = model, effect=effect)
se4.911 <- coeftest(m4.911, vcov.=cclust)[,2] 
W4.911 <- linearHypothesis(m4.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.d.911 <- plm(f2d, subset(Tscs, p911==1), index = c("COW","year"),  model = model, effect=effect)
se2.d.911 <- coeftest(m2.d.911, vcov.=cclust)[,2] 
W2.d.911 <- linearHypothesis(m2.d.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.d.911 <- plm(f2d, data=subset(Tscs, oecd==0 & p911==1), index = c("COW","year"),  model = model, effect=effect)
se3.d.911 <- coeftest(m3.d.911, vcov.=cclust)[,2] 
W3.d.911 <- linearHypothesis(m3.d.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.d.911 <- plm(f2d, data=subset(Tscs, oecd==1 & p911==1), index = c("COW","year"),  model = model, effect=effect)
se4.d.911 <- coeftest(m4.d.911, vcov.=cclust)[,2] 
W4.d.911 <- linearHypothesis(m4.d.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.r.911 <- plm(f2r, subset(Tscs, p911==1), index = c("COW","year"),  model = model, effect=effect)
se2.r.911 <- coeftest(m2.r.911, vcov.=cclust)[,2] 
W2.r.911 <- linearHypothesis(m2.r.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.r.911 <- plm(f2r, data=subset(Tscs, oecd==0 & p911==1), index = c("COW","year"),  model = model, effect=effect)
se3.r.911 <- coeftest(m3.r.911, vcov.=cclust)[,2] 
W3.r.911 <- linearHypothesis(m3.r.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.r.911 <- plm(f2r, data=subset(Tscs, oecd==1 & p911==1), index = c("COW","year"),  model = model, effect=effect)
se4.r.911 <- coeftest(m4.r.911, vcov.=cclust)[,2] 
W4.r.911 <- linearHypothesis(m4.r.911,  c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


Wald911 <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.911, W3.911, W4.911, W2.d.911, W3.d.911, W4.d.911, W2.r.911, W3.r.911, W4.r.911), digits = 3) )), 
                c("Country fixed effects", rep("\\text{Yes}", 9)),
                c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.911,m3.911,m4.911, m2.d.911,m3.d.911,m4.d.911 ,m2.r.911,m3.r.911,m4.r.911,
          se=list(se2.911,se3.911,se4.911, se2.d.911, se3.d.911, se4.d.911, se2.r.911,se3.r.911,se4.r.911),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=Wald911,
          covariate.labels = covlab[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Post-9/11 sample",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.911.tex"
) 


### Table A10: Muslim vs. Non-Muslim refugees
f2.mus <- update(f0, lnattack ~ lnref_mus.l1 + lnref_nmus.l1 + .)
f2.d.mus <- update(f0, lnattack_dom ~ lnref_mus.l1 + lnref_nmus.l1 + .)
f2.r.mus <- update(f0, lnattack_isref ~ lnref_mus.l1 + lnref_nmus.l1 + .)

m2.mus <- plm(f2.mus, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.mus <- coeftest(m2.mus, vcov.=cclust)[,2] 
W2.mus <- linearHypothesis(m2.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.mus <- plm(f2.mus, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.mus <- coeftest(m3.mus, vcov.=cclust)[,2] 
W3.mus <- linearHypothesis(m3.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.mus <- plm(f2.mus, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.mus <- coeftest(m4.mus, vcov.=cclust)[,2] 
W4.mus <- linearHypothesis(m4.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.d.mus <- plm(f2.d.mus, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.d.mus <- coeftest(m2.d.mus, vcov.=cclust)[,2] 
W2.d.mus <- linearHypothesis(m2.d.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.d.mus <- plm(f2.d.mus, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.d.mus <- coeftest(m3.d.mus, vcov.=cclust)[,2] 
W3.d.mus <- linearHypothesis(m3.d.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.d.mus <- plm(f2.d.mus, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.d.mus <- coeftest(m4.d.mus, vcov.=cclust)[,2] 
W4.d.mus <- linearHypothesis(m4.d.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.r.mus <- plm(f2.r.mus, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.r.mus <- coeftest(m2.r.mus, vcov.=cclust)[,2] 
W2.r.mus <- linearHypothesis(m2.r.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.r.mus <- plm(f2.r.mus, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.r.mus <- coeftest(m3.r.mus, vcov.=cclust)[,2] 
W3.r.mus <- linearHypothesis(m3.r.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.r.mus <- plm(f2.r.mus, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.r.mus <- coeftest(m4.r.mus, vcov.=cclust)[,2] 
W4.r.mus <- linearHypothesis(m4.r.mus, c("lnref_mus.l1=lnref_nmus.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

covlab.mus <- covlab
covlab.mus[2:3] <- c("ln Muslim refugees ($\\beta_{1}$)", "ln non-Muslim refugees ($\\beta_{2}$)")

WaldMUSLIM <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.mus, W3.mus, W4.mus, W2.d.mus, W3.d.mus, W4.d.mus, W2.r.mus, W3.r.mus, W4.r.mus), digits = 3) )), 
                   c("Country fixed effects", rep("\\text{Yes}", 9)),
                   c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.mus,m3.mus,m4.mus, m2.d.mus,m3.d.mus,m4.d.mus ,m2.r.mus,m3.r.mus,m4.r.mus,
          se=list(se2.mus,se3.mus,se4.mus, se2.d.mus, se3.d.mus, se4.d.mus, se2.r.mus,se3.r.mus,se4.r.mus),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldMUSLIM,
          covariate.labels = covlab.mus[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Muslim vs. Non-muslim refugees",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.MUSLIM.tex"
) 


### Table A11: Terrorism vs. civil war in sending state

f2.tntcw <- update(f0, lnattack ~ lnref_tnt_cw.l1 + lnref_tnt_ncw.l1 + lnref_ntnt_cw.l1 + lnref_ntnt_ncw.l1 + .)
f2.d.tntcw <- update(f0, lnattack_dom ~ lnref_tnt_cw.l1 + lnref_tnt_ncw.l1 + lnref_ntnt_cw.l1 + lnref_ntnt_ncw.l1 + .)
f2.r.tntcw <- update(f0, lnattack_isref ~ lnref_tnt_cw.l1 + lnref_tnt_ncw.l1 + lnref_ntnt_cw.l1 + lnref_ntnt_ncw.l1 + .)

covlab.tntcw <- c("ln Refugees",
                  "ln Refugees w/ CW  w/ TNT ($\\beta_{1}$)", "ln Refugees w/ CW  w/o TNT ($\\beta_{2}$)", 
                  "ln Refugees w/o CW  w/ TNT ($\\beta_{3}$)", "ln Refugees w/o CW  w/o TNT ($\\beta_{4}$)", 
                  "ln GDP p.c.", "ln Population", 
                  "Polity2", "Polity2$^{2}$", 
                  "Civil war", 
                  "Interstate war", 
                  "Intervention", "Neighboring civil war"
                  , "Neighboring TNT"
)
k.tntcw <- length(covlab.tntcw)

m2.tntcw <- plm(f2.tntcw, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.tntcw <- coeftest(m2.tntcw, vcov.=cclust)[,2] 
W2.tntcw <- linearHypothesis(m2.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.tntcw <- plm(f2.tntcw, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.tntcw <- coeftest(m3.tntcw, vcov.=cclust)[,2] 
W3.tntcw <- linearHypothesis(m3.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.tntcw <- plm(f2.tntcw, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.tntcw <- coeftest(m4.tntcw, vcov.=cclust)[,2] 
W4.tntcw <- linearHypothesis(m4.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.d.tntcw <- plm(f2.d.tntcw, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.d.tntcw <- coeftest(m2.d.tntcw, vcov.=cclust)[,2] 
W2.d.tntcw <- linearHypothesis(m2.d.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.d.tntcw <- plm(f2.d.tntcw, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.d.tntcw <- coeftest(m3.d.tntcw, vcov.=cclust)[,2] 
W3.d.tntcw <- linearHypothesis(m3.d.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.d.tntcw <- plm(f2.d.tntcw, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.d.tntcw <- coeftest(m4.d.tntcw, vcov.=cclust)[,2] 
W4.d.tntcw <- linearHypothesis(m4.d.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.r.tntcw <- plm(f2.r.tntcw, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.r.tntcw <- coeftest(m2.r.tntcw, vcov.=cclust)[,2] 
W2.r.tntcw <- linearHypothesis(m2.r.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.r.tntcw <- plm(f2.r.tntcw, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.r.tntcw <- coeftest(m3.r.tntcw, vcov.=cclust)[,2] 
W3.r.tntcw <- linearHypothesis(m3.r.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.r.tntcw <- plm(f2.r.tntcw, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.r.tntcw <- coeftest(m4.r.tntcw, vcov.=cclust)[,2] 
W4.r.tntcw <- linearHypothesis(m4.r.tntcw, c("lnref_tnt_cw.l1=lnref_ntnt_cw.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


Waldtntcw <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.tntcw, W3.tntcw, W4.tntcw, W2.d.tntcw, W3.d.tntcw, W4.d.tntcw, W2.r.tntcw, W3.r.tntcw, W4.r.tntcw), digits = 3) )), 
                  c("Country fixed effects", rep("\\text{Yes}", 9)),
                  c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.tntcw,m3.tntcw,m4.tntcw, m2.d.tntcw,m3.d.tntcw,m4.d.tntcw ,m2.r.tntcw,m3.r.tntcw,m4.r.tntcw,
          se=list(se2.tntcw,se3.tntcw,se4.tntcw, se2.d.tntcw, se3.d.tntcw, se4.d.tntcw, se2.r.tntcw,se3.r.tntcw,se4.r.tntcw),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=Waldtntcw,
          covariate.labels = covlab.tntcw[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Terrorism vs. civil war in sending state",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.tntcw.tex"
) 


### TAble A12: lagged dependent variables
f2.ldv <- update(f0, lnattack ~ lnref_tnt.l1 + lnref_ntnt.l1 + lnattack.l1 +  .)
f2.d.ldv <- update(f0, lnattack_dom ~ lnref_tnt.l1 + lnref_ntnt.l1 + lnattack_dom.l1 + .)
f2.r.ldv <- update(f0, lnattack_isref ~ lnref_tnt.l1 + lnref_ntnt.l1 + lnattack_isref.l1 + .)


covlab.ldv <- c("ln Refugees",
                "ln Refugees w/ TNT ($\\beta_{1}$)", "ln Refugees w/o TNT ($\\beta_{2}$)", 
                "ln Transnational attacks, $t-1$", 
                "ln Domestic attacks, $t-1$", 
                "ln Attacks against refugees, $t-1$", 
                "ln GDP p.c.", "ln Population", 
                "Polity2", "Polity2$^{2}$", 
                "Civil war", 
                "Interstate war", 
                "Intervention", "Neighboring civil war"
                , "Neighboring TNT"
)
k.ldv <- length(covlab.ldv)


m2.ldv <- plm(f2.ldv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.ldv <- coeftest(m2.ldv, vcov.=cclust)[,2] 
W2.ldv <- linearHypothesis(m2.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.ldv <- plm(f2.ldv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.ldv <- coeftest(m3.ldv, vcov.=cclust)[,2] 
W3.ldv <- linearHypothesis(m3.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.ldv <- plm(f2.ldv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.ldv <- coeftest(m4.ldv, vcov.=cclust)[,2] 
W4.ldv <- linearHypothesis(m4.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.d.ldv <- plm(f2.d.ldv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.d.ldv <- coeftest(m2.d.ldv, vcov.=cclust)[,2] 
W2.d.ldv <- linearHypothesis(m2.d.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.d.ldv <- plm(f2.d.ldv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.d.ldv <- coeftest(m3.d.ldv, vcov.=cclust)[,2] 
W3.d.ldv <- linearHypothesis(m3.d.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.d.ldv <- plm(f2.d.ldv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.d.ldv <- coeftest(m4.d.ldv, vcov.=cclust)[,2] 
W4.d.ldv <- linearHypothesis(m4.d.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.r.ldv <- plm(f2.r.ldv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.r.ldv <- coeftest(m2.r.ldv, vcov.=cclust)[,2] 
W2.r.ldv <- linearHypothesis(m2.r.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.r.ldv <- plm(f2.r.ldv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.r.ldv <- coeftest(m3.r.ldv, vcov.=cclust)[,2] 
W3.r.ldv <- linearHypothesis(m3.r.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.r.ldv <- plm(f2.r.ldv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.r.ldv <- coeftest(m4.r.ldv, vcov.=cclust)[,2] 
W4.r.ldv <- linearHypothesis(m4.r.ldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


Waldldv <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.ldv, W3.ldv, W4.ldv, W2.d.ldv, W3.d.ldv, W4.d.ldv, W2.r.ldv, W3.r.ldv, W4.r.ldv), digits = 3) )), 
                c("Country fixed effects", rep("\\text{Yes}", 9)),
                c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.ldv,m3.ldv,m4.ldv, m2.d.ldv,m3.d.ldv,m4.d.ldv ,m2.r.ldv,m3.r.ldv,m4.r.ldv,
          se=list(se2.ldv,se3.ldv,se4.ldv, se2.d.ldv, se3.d.ldv, se4.d.ldv, se2.r.ldv,se3.r.ldv,se4.r.ldv),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=Waldldv,
          covariate.labels = covlab.ldv[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Models with lagged dependent variables",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.ldv.tex"
) 


### Table A13: comprehensive lagged dependent variables
f2.cldv <- update(f0, lnattack ~ lnref_tnt.l1 + lnref_ntnt.l1 + lnattack.l1 + lnattack_dom.l1 + lnattack_isref.l1 +  .)
f2.d.cldv <- update(f0, lnattack_dom ~ lnref_tnt.l1 + lnref_ntnt.l1 + lnattack.l1 + lnattack_dom.l1 + lnattack_isref.l1 + .)
f2.r.cldv <- update(f0, lnattack_isref ~ lnref_tnt.l1 + lnref_ntnt.l1 + lnattack.l1 + lnattack_dom.l1 + lnattack_isref.l1 + .)


covlab.cldv <- c("ln Refugees",
                 "ln Refugees w/ TNT ($\\beta_{1}$)", "ln Refugees w/o TNT ($\\beta_{2}$)", 
                 "ln Transnational attacks, $t-1$", 
                 "ln Domestic attacks, $t-1$", 
                 "ln Attacks against refugees, $t-1$", 
                 "ln GDP p.c.", "ln Population", 
                 "Polity2", "Polity2$^{2}$", 
                 "Civil war", 
                 "Interstate war", 
                 "Intervention", "Neighboring civil war"
                 , "Neighboring TNT"
)
k.cldv <- length(covlab.cldv)

m2.cldv <- plm(f2.cldv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.cldv <- coeftest(m2.cldv, vcov.=cclust)[,2] 
W2.cldv <- linearHypothesis(m2.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.cldv <- plm(f2.cldv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.cldv <- coeftest(m3.cldv, vcov.=cclust)[,2] 
W3.cldv <- linearHypothesis(m3.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.cldv <- plm(f2.cldv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.cldv <- coeftest(m4.cldv, vcov.=cclust)[,2] 
W4.cldv <- linearHypothesis(m4.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.d.cldv <- plm(f2.d.cldv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.d.cldv <- coeftest(m2.d.cldv, vcov.=cclust)[,2] 
W2.d.cldv <- linearHypothesis(m2.d.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.d.cldv <- plm(f2.d.cldv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.d.cldv <- coeftest(m3.d.cldv, vcov.=cclust)[,2] 
W3.d.cldv <- linearHypothesis(m3.d.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.d.cldv <- plm(f2.d.cldv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.d.cldv <- coeftest(m4.d.cldv, vcov.=cclust)[,2] 
W4.d.cldv <- linearHypothesis(m4.d.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m2.r.cldv <- plm(f2.r.cldv, data=Tscs, index = c("COW","year"),  model = model, effect=effect)
se2.r.cldv <- coeftest(m2.r.cldv, vcov.=cclust)[,2] 
W2.r.cldv <- linearHypothesis(m2.r.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m3.r.cldv <- plm(f2.r.cldv, data=subset(Tscs, oecd==0), index = c("COW","year"),  model = model, effect=effect)
se3.r.cldv <- coeftest(m3.r.cldv, vcov.=cclust)[,2] 
W3.r.cldv <- linearHypothesis(m3.r.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]

m4.r.cldv <- plm(f2.r.cldv, data=subset(Tscs, oecd==1), index = c("COW","year"),  model = model, effect=effect)
se4.r.cldv <- coeftest(m4.r.cldv, vcov.=cclust)[,2] 
W4.r.cldv <- linearHypothesis(m4.r.cldv, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=cclust)$"Pr(>Chisq)"[2]


Waldcldv <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2.cldv, W3.cldv, W4.cldv, W2.d.cldv, W3.d.cldv, W4.d.cldv, W2.r.cldv, W3.r.cldv, W4.r.cldv), digits = 3) )), 
                 c("Country fixed effects", rep("\\text{Yes}", 9)),
                 c("Year fixed effects",    rep("\\text{Yes}", 9)) )
stargazer(m2.cldv,m3.cldv,m4.cldv, m2.d.cldv,m3.d.cldv,m4.d.cldv ,m2.r.cldv,m3.r.cldv,m4.r.cldv,
          se=list(se2.cldv,se3.cldv,se4.cldv, se2.d.cldv, se3.d.cldv, se4.d.cldv, se2.r.cldv,se3.r.cldv,se4.r.cldv),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=Waldcldv,
          covariate.labels = covlab.cldv[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Models with lagged attacks",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.cldv.tex"
) 


### Table A14: dyadic analyses
f1.dy <- Formula( log(attack+1) 
                  ~ log(refugees_sum+1)*attacktnhomed
                  + lngdp.1 + lnpop.1 +pol2.1 + pol2sq.1 + cw.1 + Wcw.d + Wtnt.d 
                  + lngdp.2 + lnpop.2 +pol2.2 + pol2sq.2 + cw.2 
                  + war12 +intervention   |  year+ dyadid  | 0 | dyadid)    

m1a.dy <- felm(f1.dy, data=subset(ddata1, oecd.1==0), na.action = na.omit)
m1b.dy <- felm(f1.dy, data=subset(ddata1, oecd.1==1), na.action = na.omit)

f3.dy <- update(f1.dy, log(attack_isref+1) ~ .)
m3a.dy <- felm(f3.dy, data=subset(ddata1, oecd.1==0), na.action = na.omit)
m3b.dy <- felm(f3.dy, data=subset(ddata1, oecd.1==1), na.action = na.omit)

k <- length(m1a.dy$beta)
clab <- c("ln Refugees$_{h}$", "TNT$_{s}$", "ln Refugees$_{h}$ $\\times$ TNT$_{s}$", 
          "ln GDP p.c.$_{h}$", "ln Population$_{h}$", "Polity2$_{h}$", "Polity2$^{2}$$_{h}$", "Civil war$_{h}$", "Neighboring civil war$_{h}$", "Neighboring TNT$_{h}$", 
          "ln GDP p.c.$_{s}$", "ln Population$_{s}$", "Polity2$_{s}$", "Polity2$^{2}$$_{s}$", "Civil war$_{s}$", "Interstate war$_{sh}$", "Intervention$_{sh}$")


stargazer(m1a.dy, m1b.dy, m3a.dy, m3b.dy,
          omit.stat=c("f",  "ser", "adj.rsq"),
          covariate.labels = clab,
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          add.lines = list(   c("Dyad fixed effects", "Yes", "Yes","Yes","Yes"),
                              c("Year fixed effects", "Yes", "Yes","Yes","Yes")),
          dep.var.labels = c("Transnational attacks",  "Attacks against refugees"),
          title = "Dyadic Models",
          column.labels = c( rep(c("non-OECD", "OECD"),2)),
          order=c(1:2, k, 3:(k-1)),
          notes = "Clustered standard errors in parentheses",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE,
          out= "table.dyadic.tex"
)


### Table A15: Quasi Poisson

sec <- function(x) vcovHC(x, cluster="COW", type="HC0")

m2 <- glm(update(f2, attack ~ . + as.factor(year)), data=Tscs, family=quasipoisson)
se2 <- coeftest(m2, vcov.=sec)[,2] 
W2 <- linearHypothesis(m2, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m3 <- glm(update(f2, attack ~ . + as.factor(year)), data=subset(Tscs, oecd==0), family=quasipoisson)
se3 <- coeftest(m3, vcov.=sec)[,2] 
W3 <- linearHypothesis(m3, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m4 <- glm(update(f2, attack ~ . + as.factor(year)), data=subset(Tscs, oecd==1), family=quasipoisson)
se4 <- coeftest(m4, vcov.=sec)[,2] 
W4 <- linearHypothesis(m4, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]


m2d <- glm(update(f2d, attack_domestic ~ . + as.factor(year)), data=Tscs, family=quasipoisson)
se2d <- coeftest(m2d, vcov.=sec)[,2]
W2d <- linearHypothesis(m2d, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m3d <- glm(update(f2d, attack_domestic ~ . + as.factor(year)), data=subset(Tscs, oecd==0), family=quasipoisson)
se3d <- suppressWarnings(coeftest(m3d, vcov.=sec)[,2])
W3d <- linearHypothesis(m3d, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m4d <- glm(update(f2d, attack_domestic ~ . + as.factor(year)), data=subset(Tscs, oecd==1), family=quasipoisson)
se4d <- suppressWarnings(coeftest(m4d, vcov.=sec)[,2])
W4d <- linearHypothesis(m4d, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m2r <- glm(update(f2r, attack_isref ~ . + as.factor(year)), data=Tscs, family=quasipoisson)
se2r <- suppressWarnings(coeftest(m2r, vcov.=sec)[,2])
W2r <- linearHypothesis(m2r, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m3r <- glm(update(f2r, attack_isref ~ . + as.factor(year)), data=subset(Tscs, oecd==0), family=quasipoisson)
se3r <- suppressWarnings(coeftest(m3r, vcov.=sec)[,2])
W3r <- linearHypothesis(m3r, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]

m4r <- glm(update(f2r, attack_isref ~ . + as.factor(year)), data=subset(Tscs, oecd==1), family=quasipoisson)
se4r <- suppressWarnings(coeftest(m4r, vcov.=sec)[,2])
W4r <- linearHypothesis(m4r, c("lnref_tnt.l1=lnref_ntnt.l1") , vcov=sec)$"Pr(>Chisq)"[2]


WaldFE <- list(c("Wald test $\\beta_{1}=\\beta_{2}$", sprintf("%.3f", round(c(W2, W3, W4, W2d, W3d, W4d, W2r, W3r, W4r), digits = 3) ))
               , 
               c("Country fixed effects", rep("\\text{No}", 9)),
               c("Year fixed effects",    rep("\\text{Yes}", 9)) 
)
stargazer(m2,m3,m4, m2d,m3d,m4d ,m2r,m3r,m4r,
          se=list(se2,se3,se4, se2d, se3d, se4d, se2r,se3r,se4r),
          omit.stat=c("f", "adj.rsq"),
          omit="y",
          add.lines=WaldFE,
          covariate.labels = covlab[-1],
          dep.var.labels = c("Transnational attacks", "Domestic attacks", "Attacks against refugees"),
          star.cutoffs = c(0.05, 0.01, 0.001),
          type="latex",
          title = "Quasi Poisson regression estimates",
          column.labels = c( rep(c("all countries", "non-OECD", "OECD"),3)),
          notes = "Clustered standard errors in parentheses",
          float.env= "sidewaystable",
          font.size="scriptsize",
          digits=3,
          digits.extra=3,
          column.sep.width = "0pt",
          align = TRUE ,
          out= "table.COUNT.tex"
) 

# Figures A2: density plots
pdf(file ="ref_dens.pdf", width=4, height=4)
plot(density(Tscs$lnref_tnt.l1), las=1, bty='n', xlab="ln Refugees", main="all countries")
lines(density(Tscs$lnref_ntnt.l1), lty=2)
legend("topright", c("Refugees w/ TNT", "Refugees w/o TNT"), lty=1:2, bty="n")
dev.off()

pdf(file ="ref_dens_oecd.pdf", width=4, height=4)
plot(density(Tscs$lnref_tnt.l1[Tscs$oecd==0]), las=1, bty='n', xlab="ln Refugees", main="non-OECD")
lines(density(Tscs$lnref_ntnt.l1[Tscs$oecd==0]), lty=2)
legend("topright", c("Refugees w/ TNT", "Refugees w/o TNT"), lty=1:2, bty="n")
dev.off()

pdf(file ="ref_dens_noecd.pdf", width=4, height=4)
plot(density(Tscs$lnref_tnt.l1[Tscs$oecd==1]), las=1, bty='n', xlab="ln Refugees", main="OECD")
lines(density(Tscs$lnref_ntnt.l1[Tscs$oecd==1]),  lty=2)
legend("topright", c("Refugees w/ TNT", "Refugees w/o TNT"), lty=1:2, bty="n")
dev.off()


### Figures A3 and A4: Settled vs. Resettled Refugees
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  require(grid)
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}


# resettlement data (from UNHCR, cleaned)
load("resettlement2018.RData")
res <- subset(resettlement2018, Year > 1969)

# collapse by ccode1
res1 <- summaryBy(Value ~ ccode1 + Year, data = res, FUN= sum)

# refugees data (from UNHCR, cleaned)
load("refugeesonly2018_dyadic.RData")

# collapse by ccode1
ref1 <- summaryBy(Value.sum ~ ccode1 + Year, data = refugeesonly2018_dyadic, FUN= sum)


### Figure A4: Individual country plots
# United Kingdom
ukref <- ref1[ref1$ccode1 == 200, ]
ukres <- res1[res1$ccode1 == 200, ]
ukref$status <- "refugees"
ukres$status <- "resettled"
names(ukref) <- c("ccode1", "Year", "Value", "status")
names(ukres) <- c("ccode1", "Year", "Value", "status")
ukres <- ukres[complete.cases(ukres), ]
ukref <- ukref[complete.cases(ukref), ]
uk <- merge(ukref, ukres, by=c("ccode1","Year", "Value", "status"), all =T)
uk <- data.frame(uk)
uk <- subset(uk, Year > 1987)
uk <- ggplot(uk, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("United Kingdom")

# United States
usref <- ref1[ref1$ccode1 == 2, ]
usres <- res1[res1$ccode1 == 2, ]
usref$status <- "refugees"
usres$status <- "resettled"
names(usref) <- c("ccode1", "Year", "Value", "status")
names(usres) <- c("ccode1", "Year", "Value", "status")
usres <- usres[complete.cases(usres), ]
usref <- usref[complete.cases(usref), ]
us <- merge(usref, usres, by=c("ccode1","Year", "Value", "status"), all =T)
us <- data.frame(us)
us <- subset(us, Year > 1987)
us<- ggplot(us, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("United States")

#France
franref <- ref1[ref1$ccode1 == 220, ]
franres <- res1[res1$ccode1 == 220, ]
franref$status <- "refugees"
franres$status <- "resettled"
names(franref) <- c("ccode1", "Year", "Value", "status")
names(franres) <- c("ccode1", "Year", "Value", "status")
franres <- franres[complete.cases(franres), ]
franref <- franref[complete.cases(franref), ]
france <- merge(franref, franres, by=c("ccode1","Year", "Value", "status"), all =T)
france <- data.frame(france)
france <- subset(france, Year > 1987)
f<- ggplot(france, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("France")

#Germany
gerref <- ref1[ref1$ccode1 == 255, ]
gerres <- res1[res1$ccode1 == 255, ]
gerref$status <- "refugees"
gerres$status <- "resettled"
names(gerref) <- c("ccode1", "Year", "Value", "status")
names(gerres) <- c("ccode1", "Year", "Value", "status")
gerres <- gerres[complete.cases(gerres), ]
gerref <- gerref[complete.cases(gerref), ]
germany <- merge(gerref, gerres, by=c("ccode1","Year", "Value", "status"), all =T)
germany <- data.frame(germany)
germany <- subset(germany, Year > 1987)
ge <- ggplot(germany, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("Germany")

#Italy
itaref <- ref1[ref1$ccode1 == 325, ]
itares <- res1[res1$ccode1 == 325, ]
itaref$status <- "refugees"
itares$status <- "resettled"
names(itaref) <- c("ccode1", "Year", "Value", "status")
names(itares) <- c("ccode1", "Year", "Value", "status")
itares <- itares[complete.cases(itares), ]
itaref <- itaref[complete.cases(itaref), ]
italy <- merge(itaref, itares, by=c("ccode1","Year", "Value", "status"), all =T)
italy <- data.frame(italy)
italy <- subset(italy, Year > 1987)
i <- ggplot(italy, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("Italy")

# Greece
greref <- ref1[ref1$ccode1 == 350, ]
greres <- res1[res1$ccode1 == 350, ]
greref$status <- "refugees"
greres$status <- "resettled"
names(greref) <- c("ccode1", "Year", "Value", "status")
names(greres) <- c("ccode1", "Year", "Value", "status")
greres <- greres[complete.cases(greres), ]
greref <- greref[complete.cases(greref), ]
greece <- merge(greref, greres, by=c("ccode1","Year", "Value", "status"), all =T)
greece <- data.frame(greece)
greece <- subset(greece, Year > 1987)
g <- ggplot(greece, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("Greece")

# Arrange in single plot
pdf("fgA4.pdf")
multiplot(uk, us, f, ge, i, g,
          cols = 2)
dev.off()


### Figures A3: All OECD plot
# Select oECD countries from refugee data
ref1$oecd <- 0
ref1$oecd[ref1$ccode1==900 & ref1$Year>=1971 |
            ref1$ccode1==305 & ref1$Year>=1961 |
            ref1$ccode1==211 & ref1$Year>=1961 |
            ref1$ccode1==20  & ref1$Year>=1961 |
            ref1$ccode1==155 & ref1$Year>=2010 |
            ref1$ccode1==316 & ref1$Year>= 1995|
            ref1$ccode1==390 & ref1$Year>=1961 |
            ref1$ccode1==366 & ref1$Year>=2010 |
            ref1$ccode1==375 & ref1$Year>=1969 |
            ref1$ccode1==220 & ref1$Year>=1961 |
            (ref1$ccode1==255 | ref1$ccode1==260) & ref1$Year>=1961 |
            ref1$ccode1==350 & ref1$Year>=1961 |
            ref1$ccode1==310 & ref1$Year>=1996 |
            ref1$ccode1==395 & ref1$Year>=1961 |
            ref1$ccode1==205 & ref1$Year>=1961 |
            ref1$ccode1==666 & ref1$Year>=2010 |
            ref1$ccode1==325 & ref1$Year>=1962 |
            ref1$ccode1==740 & ref1$Year>=1964 |
            ref1$ccode1==730 & ref1$Year>=1996 |
            ref1$ccode1==367 & ref1$Year>=2016 |
            ref1$ccode1==212 & ref1$Year>=1961 |
            ref1$ccode1==70  & ref1$Year>=1994 |
            ref1$ccode1==210 & ref1$Year>=1961 |
            ref1$ccode1==920 & ref1$Year>=1973 |
            ref1$ccode1==385 & ref1$Year>=1961 |
            ref1$ccode1==290 & ref1$Year>=1996 |
            ref1$ccode1==235 & ref1$Year>=1961 |
            ref1$ccode1==317 & ref1$Year>=2000 |
            ref1$ccode1==349 & ref1$Year>=2010 |
            ref1$ccode1==230 & ref1$Year>=1961 |
            ref1$ccode1==380 & ref1$Year>=1961 |
            ref1$ccode1==225 & ref1$Year>=1961 |
            ref1$ccode1==640 & ref1$Year>=1961 |
            ref1$ccode1==200 & ref1$Year>=1961 |
            ref1$ccode1==2 & ref1$Year>=1961 ] <- 1
ref2 <- subset(ref1, oecd == 1)

# aggregate by year
ref2 <- summaryBy(Value.sum.sum ~ Year, data = ref2, FUN= sum)
ref2$status <- "refugees"
names(ref2) <- c("Year", "Value", "status")



# Select oECD countries from resettlement data
res1$oecd <- 0
res1$oecd[  res1$ccode1==900 & res1$Year>=1971 |
              res1$ccode1==305 & res1$Year>=1961 |
              res1$ccode1==211 & res1$Year>=1961 |
              res1$ccode1==20  & res1$Year>=1961 |
              res1$ccode1==155 & res1$Year>=2010 |
              res1$ccode1==316 & res1$Year>= 1995|
              res1$ccode1==390 & res1$Year>=1961 |
              res1$ccode1==366 & res1$Year>=2010 |
              res1$ccode1==375 & res1$Year>=1969 |
              res1$ccode1==220 & res1$Year>=1961 |
              (res1$ccode1==255 | res1$ccode1==260) & res1$Year>=1961 |
              res1$ccode1==350 & res1$Year>=1961 |
              res1$ccode1==310 & res1$Year>=1996 |
              res1$ccode1==395 & res1$Year>=1961 |
              res1$ccode1==205 & res1$Year>=1961 |
              res1$ccode1==666 & res1$Year>=2010 |
              res1$ccode1==325 & res1$Year>=1962 |
              res1$ccode1==740 & res1$Year>=1964 |
              res1$ccode1==730 & res1$Year>=1996 |
              res1$ccode1==367 & res1$Year>=2016 |
              res1$ccode1==212 & res1$Year>=1961 |
              res1$ccode1==70  & res1$Year>=1994 |
              res1$ccode1==210 & res1$Year>=1961 |
              res1$ccode1==920 & res1$Year>=1973 |
              res1$ccode1==385 & res1$Year>=1961 |
              res1$ccode1==290 & res1$Year>=1996 |
              res1$ccode1==235 & res1$Year>=1961 |
              res1$ccode1==317 & res1$Year>=2000 |
              res1$ccode1==349 & res1$Year>=2010 |
              res1$ccode1==230 & res1$Year>=1961 |
              res1$ccode1==380 & res1$Year>=1961 |
              res1$ccode1==225 & res1$Year>=1961 |
              res1$ccode1==640 & res1$Year>=1961 |
              res1$ccode1==200 & res1$Year>=1961 |
              res1$ccode1==2 & res1$Year>=1961 ] <- 1
res2 <- subset(res1, oecd == 1)
res2 <- summaryBy(Value.sum ~ Year, data = res2, FUN= sum)
res2$status <- "resettled"
names(res2) <- c("Year", "Value", "status")
all <- merge(ref2, res2, by =c("Year", "Value", "status"), all =T)
all <- subset(all, Year > 1987)
all <- all[complete.cases(all), ]

pdf("fgA3.pdf")
ggplot(all, aes(fill=status, y=Value, x=Year)) + 
  geom_bar(position="stack", stat="identity") +
  ggtitle("Refugees and resettled refugees in OECD countries 1988-2018")
dev.off()
